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Abstract. The interest to retrial queueing systems is due to their application 
to telephone systems. The paper studies multiserver retrial queueing systems 
with n servers. Arrival process is a quite general point process. An arriving 
customer occupies one of free servers. If upon arrival all servers are busy, 
then the customer waits for his service in orbit, and after random time retries 
more and more to occupy a server. The orbit has one waiting space only, and 
arriving customer, who finds all servers busy and the waiting space occupied, 
losses from the system. Time intervals between possible retrials are assumed 
to have arbitrary distribution (the retrial scheme is exactly explained in the 
paper). The paper provides analysis of this system. Specifically the paper 
studies optimal number of servers to decrease the loss proportion to a given 
value. The representation obtained for loss proportion enables us to solve the 
problem numerically. The algorithm for numerical solution includes effective 
simulation, which meets the challenge of rare events problem in simulation. 
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1. Introduction 

The paper studies multiserver queueing systems with retrials and losses. Retrial 
queueing systems are an object of investigation in numerous papers (e.g. see reviews 
Artalejo [2], Artalejo and Falin [3], Falin 201 )■ The interest to these systems is due 
to their application to telephone systems. 

During the last years there has been an especial interest to queueing systems with 
retrial and losses. The different queueing systems of this type has been studied in 
0, 23 > 23 > and other papers. 

The most of known studies in this direction are devoted to analytical solution 
of one or other problem. Bocharov et al jH] studied multiclass single-server M/G/l 
queueing system with finite buffer and limited number of places in orbit. They de- 
duced relationship between steady state distributions of this system and the similar 
system with only one class of customers. Kovalenko 22] studied the loss probability 
in M/G/n retrials queueing system with two customer classes. Mandelbaum et al 
|17| studied queueing system with retrials and losses with time-dependent param- 
eters in framework of Markovian service network (see . They provide fluid and 
diffusion approximations to obtain numerical characteristics for queue-length and 
virtual waiting time processes. Atencia and Phong 0] and Bocharov, Phong and 
Atencia [7] also studied characteristics of queueing systems with retrial and losses 
with the aid of analytic methods. 

However explicit representations can be deduced for only restricted classes of 
retrial queueing systems, and application of these results for real systems is very 
difficult. 

Recently Abramov pQ studied the main multiserver retrial queueing system, 
where input stream was a quite general point process (see pQ for more details). 
The martingale approach, that has been used there, can be adapted to systems 
with retrials and losses as well, and the present paper is just devoted to analysis of 
these systems. The exact description of the system is as follows. 
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• The arrival process Alt) is a point process (the assumptions on this process 
are given later). 

• There are n servers, and an arriving customer occupies one of free servers. 

• If upon arrival all servers are busy, but the secondary queue, orbit, having only 
one space is free, then the customer occupies the orbit and retries more and more 
to occupy a server. 

• A customer, who upon arrival finds all servers busy and the orbit occupied, 
loses from the system without delay. 

• A service time of each customer is exponentially distributed random variable 
with parameter \x. 

• A time between retrials is arbitrary distributed. It is generated by point process 
D(t) (the details on this assumption are given later). The processes D{t) and A{t) 
are assumed to be disjoint (that is the probability of simultaneous jump of these 
processes is equal to 0). 

Notice, the significant difference between the main multiserver model of PP and 
the model considered in the present paper is also that retrial times in the present 
model are generally distributed, while in pQ retrial times of each customer in orbit 
are exponentially distributed. 

There is also difference in the assumptions. In the case of the main multiserver 
retrial queueing systems considered in pQ the point process A(t) must satisfy the 
condition 



leading to the stability of the system. 

In the case of the system considered in the present paper the number of customers 
in the systems is always bounded, and the stability condition is not longer 

required, and the following more weaker conditions of convergence are used 



(1.1) 




(1.2) 




and 



(1.3) 
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i.e. A, A, fi and n are not allied to one another. 

Conditions i|I.2() and Ijl.^jl are technical conditions, that are then used by the 
methods of our analysis. (We often use the Lebesgue theorem on dominated con- 
vergence.) 

The main aim of this paper is the numerical solution of a concrete problem 
associated with specific performance measure. The main question that we want to 
answer in this paper is under what number of servers the loss proportion will be 
less than a given small value? The system that will be studied in this paper is 
described by quite general point processes, and explicit analytic solution for this 
system in general case is impossible. At the same time the direct simulation of 
this system in order to answer the aforementioned question is not available either, 
since as the loss probability is chosen very small, we should deal with simulating 
rare event problem. The simulating rare events problem is a well-known problem 
having notable attention in the literature (e.g. EH] an d many others). Then 

the main result of this paper is an answer the question how to simulate this system? 

The main results of the present paper are represented by relation l|3.8|l and The- 
orem fni These results enable us to effectively simulate models to study their most 
significant performance characteristic, the loss proportion. By effective simulation 
we mean such an approach that make the challenge of the rare events problem in 
simulation. 

In our case, we discuss three relations for loss proportions: l|3.8[) . (|3.9|l and l|3.I0|l . 
Relations and (|3~TT)|) are not effective. For example, l|3.I OJI requires to sim- 

ulate the events {Qi(s— ) = n) and {Q2{s~ ) = 1} directly. These events become 
rare as n increases indefinitely. In contrast relation (|3.8[l is effective. Its expression 
contains lim^oo i _1 E L Q(s)ds, which according to Theorem l4.1l can be easily ap- 
proximated, so that the probabilities of rare events presented in calculation of this 
expression can be removed. Specifically, there can be used a change of the Poisson 
measure in (|3.2() for calculation of (|3.8|l . i.e. replace one Poisson measure to another 
one and use the Radon-Nikodim derivative. In other words, the Poisson process 
under which the required characteristic of the system has a small probability, is re- 
placed by other Poisson process under which the usual simulation becomes available, 
and then the required small probability is recalculated by using the Radon-Nikodim 
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derivative (see e.g. (22J or QH|)- Recall that increments of a Poisson process are 
exponentially distributed, and for two exponentially distributed random variables 
with parameters i?i and $2, the Radon-Nikodim derivative is j lP ' ?1 = ^ie^ 2_1?1 ^ x , 
where P# fc (x), k = 1,2, are corresponding exponential probability distributions. 

Relations and l|3.1()Jl are derivative from (12.2ft and (|2.3[1 correspondingly. 

The presence of these equations is illustrative in order to show the difference 
between direct simulation and simulation based on (|3.8|) and Theorem 14. II 

The paper is organized as follows. In Section [21 we discuss basic equations giv- 
ing us then three aforementioned relations (|3.8[) . (|3.9fl and (|3.1U|) for loss propor- 
tions. In Section |3| we use semimartingale decomposition and Lenglart-Rebollcdo 
inequality to derive then crucial relation (j3.8(l . In Section 01 we prove Theorem 14. II 
and the corollary from this theorem related to particular case where the point pro- 
cesses A(t) and D(t) are Poisson. Specifically, Theorem l4.1l establishes relationships 
for the proportions lim^oo t^ 1 f Q ¥{Qi(s) — i, Q2(s) = i}ds (i = 0, 1, . . . , n — 1; 
j = 0, 1). Sectionals devoted to analysis of the performance measure based on loss 
proportion. There is formulated the problem of optimization of the proposed per- 
formance measure and described the algorithms for its solution. Section provides 
two numerical examples. The concluding remarks are in Section 



2. Basic equations 

All point processes considered in this paper are assumed to be right-continuous 
having the left-side limits. 

The following notation is used. The arrival process is denoted A(t) . The queueing 
process (the number of occupied servers process) is denoted Q\(t). The orbit state 
is denoted Q2(t)- It can take only the values and 1 correspondingly to the cases 
of empty and occupied orbit space. The loss process is denoted Q^it). Specifically, 
Qs{t) is the cumulated number of losses up to time t. The moments of retrials are 
governed by point process D(t), that is, the retrial moments are t\, t%, ■ ■ ■ and D(t) 
= min{A: : tu < i}. The sequence of independent Poisson processes all with the 
same rate \i is denoted 7rj(£). 
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Then we have the following equations denotes the indicator of event A): 

(2.1) Qi(t) + Q 2 (t) + Q 3 (t) = A(t)- / Vl{Qi(s-)>Z}d7n(s), 

Jo 1=1 

(2.2) Q2(t)+Q 3 (t) = / HQi(s-) = n}dA(s) 

Jo 

I{Qi(s-) + n}I{Q 2 (s-) = l}dD(s), 

(2.3) Q 3 (t) = f I{Q 1 (s-) = n}I{Q 2 {s-) = l}dA(s). 

Jo 

Clearly, that the meaning of the left-hand side of relation (|2.1|) is the total number 
of customers in the system including customers in service, one in orbit (in the 
case of occupied place) as well as all losses up to time t. Notice, that the right- 
hand side of (|2.1|) is as the right-hand side of (2.1) in Q]. The first term of the 
right-hand side of i|2.2|l means the total number of arrivals to busy system where 
all servers are occupied. The second term of the right-hand side of (|2.2(l is the 
subtracted number of successfully retrial customers up to time t from the total 
number of arrivals to busy system during the same time t. Successful retrials occur 
in the cases where immediately before retrial instants there is at least one busy 
server. Notice, that relation 1)2. 2JI is structured similarly to corresponding relation 
(2.2) of 1 . Relation 1)2. 3)1 characterizes the number of losses. Losses occur in the 
cases where immediately before arrival all servers are busy and the place in orbit 
is occupied. This is indicated by the right-hand side of (|2.3|) . 

3. Semimartingale decompositions and normalization 

In this section we discuss the normalized processes by dividing these processes 
by t. As t increases to infinity, the behaviour of the process Qz{t) divided by t is 
the significant performance characteristics characterizing losses per time unit. 

For the purpose of study normalized characteristics we need in semimartin- 
gale decomposition of the processes (e.g. Liptser and Shiryayev Jacod and 
Shiryayev [TT]V 

The processes A(t), D{t) and 717(f), I — 1, 2, . . . , m all are assumed to be given on 
the common filtered probability space (fl, J 7 , F = (J 7 t)t>0: !")■ All these processes 
are semimartingales. The semimartingale decomposition for arbitrary point process 
X(t) is written X{t) = X(t) + Mx(t), where X(t) is a compensator and Mx(t) 
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is a local square integrable martingale. For example: A(t)=A(t)+MA{t). In some 
cases for local square integrable martingales we use the notation M with indexes 
such as Mij t k(t). Such types of notation will be especially explained. 
Observing first l|2.1|l . rewrite it 

(3.1) Qi(t) + Q 2 {t) + Q 3 (t) = A(t)-C(t), 

(3.2) C(t) = / Vl{Qi(s-)>Z}d7r / ( S ). 

The semimartingale C(t) admits Doob-Meyer decomposition 

(3.3) C(t) = C(t)+M c (t). 

In turn, the compensator C(t) admits the representation 

(3.4) C{t) =n f Qi(«)da 

Jo 

(for details see Dellacherie [5], Liptser and Shiryayev ^3], ^5], Theorem 1.6.1). 
Therefore, in view of (|3.2|l . 1(3.3(1 and 13.4(1 . one can rewrite 1(3. l|l as follows: 

(3.5) Qt{t) + Q 2 (t) + Q 3 (t) = A(t) -fx [ Qx(s)ds - M c {t). 

Jo 

Let us pass to normalized processes. For arbitrary process X(t), for i > its 
normalization is denoted by small letter, so 

««> = 

Therefore, passing to normalized processes in ((3.511 can we written as 

(3.6) qi (t) + q 2 (t) + q 3 (t) = a(t) - ^ Q 1 (s)ds - m c (t). 

1 Jo 

Assume now that t increases to infinity. Then, because of stochastic boundedness 
of Qi(t) and Q 2 (t) the terms qi(t) and q 2 (t) vanish with probability 1, and the 
left-hand side of ((3.6(1 in limit looks P" Hindoo q^ (t). 
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Next, according to 1|1.2|) . as t increases to infinity, a(t) converges with probabil- 
ity 1 to A. The term mc(t) vanishes in probability, since according to Lenglart- 
Rebolledo inequality we have: 

smc(s) 



\m c (t)\ > 6} < P{ sup 

^0<s<t 



> S 



(3.7) 



sup \M c (s)\ > St 

0<s<t 

and because of arbitrariness of e in l|3.7|l the statement is right. 

Therefore, applying the Lebesgue theorem on dominated convergence 

P- lim q 3 (t) = A - /i (V Um - [ Q^ds] 

t—><X \ t^OO tin J 

(3.8) V ° J 

= A - n y lim jEj Qi(s)ds^ 

The meaning of the term of the right-hand side of 1|3.8|) . taken in brackets, is 
the expected number of occupied servers during the long-run history, or expected 
stationary number of occupied servers. 

Relations <|2.2[) and (|2.3() provide additional information on the behaviour of the 
process q 3 (t). Specifically, we have 

1 f* 
im q 3 (t) = lim — E / 



(3.9) 



and 



lim q 3 (t) = lim -E / I{Qi(s-) = n}dA(s) 

t — >-00 t— >00 t Jq 

- lim jE / I{Qi(s-) ^ n}I{Q 2 {s-) = l}dD(s) 

t^oo t Jn 



1 f l 

(3.10) P" lim q 3 (t) = lim -E / I{Qi(s-) = n}I{Q 2 (s-) = l}d A(s). 

t — >CG t — >00 t Jq 

In the particular case where the process A(t) is Poisson from l|3.10|l we obtain: 
P" lim q 3 (t) = A lim -E / I{Qi(s) = n}I{Q 2 (s) = l}ds 

t^oo t— too tin 

= A lim - / P{Qi(s) = n,Q 2 (s) = l}ds, 

t^OO t Jn 
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and from l|3.9|l we obtain 



(3.12) 



i r* 

■ lim q 3 (t) =A lim - / P{Qi(s) = n}ds 

t—>OC t-^-OC t Jq 



lim \e f I{Qi(«-) 7^ n}I{Q 2 (s-) = l}d£>(a). 



In the case where both A(t) and D{t) are Poisson, qs(t) converges to its limit with 
probability 1 as t increases to infinity, because both mi(t) and mo{t) vanish with 
probability 1. Specifically, in this case we obtain: 



lim q 3 (t) = A lim - / V{Qi(s) = n}ds 

t — >OQ t — >OQ f Jq 



(3.13) 



- A lim \ / P{Qi(«) + n, Q 2 (s) = l}ds) 
1 { lim q 3 (t) = A lim P{Qi(t) = n} 

Li — > oc £— >-oc 

-A lim P{Qi(t)^n,Q 2 (*) = 1}} 



There is use of the fact of existence of the limiting stationary probabilities in (|3.13|) 
as t — > oo, as well as the equalities 

i r* 

lim P{Qi(t) =n} = lim - / P{Qi(s) = n}ds, 

£^■00 t—>-OG t Jq 

1 /"* 

lim P{Qi(i) ^ n,Q 3 (t) = 1} = lim - / P{Q x (s) ^ n,Q 2 («) = l}da. 

t — 'OO t — >00 t Jq 

4. Limiting frequencies and distributions 

The results obtained in the previous section do not provide completed infor- 
mation about losses in this system. We have no information about the quantity 
lim^oo t~ x J"' Qi(s)ds entering to the right-hand side of JITS}. The information 
about the right-hand sides of l|3. 9fl and 13.10|) is unknown as well. 

Therefore in this section we derive equations for the following frequencies: 

1 I* 

lim - / P{Q 1 (s)=i,Q 2 (s)=j}ds, 
(4.1) t Jo 

i = 0, l,...,n; j = 0,1. 
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For this purpose we introduce the processes: 

l i j(t) = l{Q 1 (t) = i n Q 2 (t)=j}, 

(4.2) 

i = 0, 1,.. . ,n; j = 0,1, 

taking the value if at least one of the indexes i or j is negative. 

The jump of a point process is denoted by adding A 1 . For example, AA(t) is a 
jump of -A(£), Ani(t) is a jump for 717 (t) and so on. 

Denote 

i 

ni(t) = £V,(t). 

Then we have the following equations: 
(4.3) I{0i(*-) + AQi(t) = i n Q 2 (i) = Q 2 (t-) = i} 

= / i _ 1;j (i-)AA(*) + / j+ i >j (t-)An i+1 (t) 

+ I id (t-)[1 - AA(t)][l - An,(t)][l - jAD(t)} 7 
i = 0, 1,.. .,n- 1; j = 0,1, 



i{Qi(t-) + AQi(t) = z n Q 2 (t-) = 1 n Q 2 (t) = 0} 

(4.4) =/<_!,! (t-)AD(t), 

i = 0, 1, . . . ,n — 1, 



(4.5) 



(4.6) 



(4.7) 



+ AQi(t) = n n Q 2 (t) = Qa(t-) = 0} 
= / n _i, (t-)AA(t) + J n ,o(*-)[l - An„(t)][l - AA(t)], 

+ AQi(t) = n n g 2 (t-) = n Q 2 (*) = 1} 
= I nfi (t-)AA(t), 

UQi(t-) + AQi(t) = n n Q 2 (i) = Q 2 (t-) = 1} 
= I n - 1A (t-)AA(t) + I nA (t-)AA(t) + /„,i(t-)[l - AA(t)][l - An„], 



lr This is the triangle sign. We hope that this sign will not be mixed with the parameter (capital 
Greek letter) A 
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I{Qi(t-) + AQi(i)=n n Q 2 (t-) = l n Q 2 {t) = 0} 

(4.8) 

= I n -i,i (t-)AD(t). 

Then, 

AIi,j(t) =I{Qi(t-) + AQi(t)=i n Q 2 (i) = Q 2 (i-)=j} 

+ I{Qi(i-)+AQi(t)=i n Q 2 (t-)^Q 2 (*)=j} 

(4.9) 

i = 0, 1, . . . ,n; j = 0, 1. 

Since 

=kj(t) 

then taking into account (|4.3I) - 14. 9|) we have the following. 
For i = 0, 1, . . . , n — 1, 

I M (i) =1^(0) + (\h-xjis-) - I tJ ( S -)]dA(s) 
Jo 

(4.10) - f I itj (s-)dILi(s) - j [ I itj (s-)dD(s) 

Jo Jo 

+ [ / i +i.j(a-)dn i+ i(s) + (l-j) / (s-)dD(a). 
Jo Jo 

For i — n we have the following pair of equations. In the case j = we have: 

Jn,o(t) =/n,o(0)+ / [J n _l, (s-) - /„, (s-)]dA(s) 







(4.ii) - / / n)0 ( s -)dn n ( s ) 







+ / I n -i,x(s-)dD(s). 
Jo 

In turn, in the case j ' = 1 we have: 

=/ «,i(0)+ / J„_i,i(«-)dA(«) 



(4.12) 





+ / I n ,o(s-)dA(s) - f / n , 1 (s-)dn n (s). 
Jo Jo 



12 VYACHESLAV M. ABRAMOV 

By semimartingale decomposition, from Ij4.1t)|l . (|4.11(l and (|4.12(l wc obtain the 
following equations. In the cases i = 0, 1, . . . , n — 1; j = 0, 1 we have: 

I id {t) =I id (0) + [ [It-iA*-) - I<j(*-)]M{8) 



-ill I Iij(s)ds + (i + / I i+1 j(s)ds 
(4.13) Jo Jo 

-j f 7 M (s-)dU(a) + (l-j) / Ii-i,i(s-)dD(s) 



with the martingale 



(4.14) Mi d (t) = - I itj (s-)dMnM+ 7 i+1 , i (s-)dM n<+1 (s). 

Jo Jo 

For i — n we have the following pair of equations. In the case j = we have: 

4,o(<) =I»,o(0) + / [/„_i, (s-) - I„, {s-)]dA(s) 
Jo 

(4 ' 15) -nil! I n>0 (s)ds+ ! I n ^ ltl (s-)dD(s) 



+ M„ t0 (t), 

with the martingale 

(4.16) M„, (t) = - f I n , {s-)dMnM- 

Jo 

In turn, in the case j = 1 we have: 



I n .i(t) =4,i(0)+ / 4- M (s-)dA( S )-np / 7«,i(*)ds 
(4.17) Jo 

+ / / n)0 («-)dA(s)+Af n ,i(t), 



with the martingale 
(4.18) M„, 1 (t) = - / / n ,i(s-)dM n „(s). 



Theorem 4.1. for i/ie limiting state frequencies of the processes Qi(t) and Q2(t) 
we have the following equations. 
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In the boundary cases i = for j = 0, 1 we have: 

i r* 

H lim - / P{Qi(«) = l,Q 2 (s)=j}ds 

(4.19) = lim 1e / I{Qi(s-) = 0, Q 2 (s-) - j}<L4(s) 

+ j lim / I{Qi(«-) = 0,Q 3 (O=i}dD(*)- 

In the cases i = 1, 2, . . . , n — 1, j = 0, 1 we /law: 
1 /■' 

i/i lim - / F{Q 1 (s)=i,Q 2 (s)=j}ds 

1 ^ 

- (i + lim - / P{Qi(s) = i + 1, Q 2 (s) = j}ds 

i Jo 

= lim ±E / [I{Qi(s-) = i - 1, Q 2 (s-) = j} 

(4.20) *- ,0 ° 1 Jo 

-I{Q 1 (s-)=i,Q 2 (s-)=j}]dA(s) 
1 /•* 

+ (l-i) lim -E / I{Q 1 (s-)=i-l ) Q 2 (s-) = l}d£> s) 
* Jo 

1 /■* 

- j lim -E / I{Qi(s-) = i,Q 2 (s-) = l}dD(s). 

In the case i — n and j — we have: 
1 f* 

n\i lim - / P{Qi(s) = n,Q 2 (s) = 0}ds 

t^oo t J 

= lim / [I{Qi(s-) = n - 1, Q 2 (a-) = 0} 

(4.21) i Jo 

-I{Qi(«-)="iQa(*-)=0}]<iA(s) 

1 r* 

+ lim -E / I{Qi(s-) = n,Q a (*-) = l}dD(«). 
In £fte case i — n and j = 1 we have: 

i r* 

nil lim - / P{Qi s = rj,Q 2 (s = l}cfe 

(4.22) = lim ]e f I{Qi(s-) = n - 1, Q 2 (s-) = 1}cM(s) 

t^oo f _y 

1 /•« 

+ lim -E / I{Qi(«-) = n, Q 2 s-) = 0}dA(s). 

Proof. The proof of this theorem easily follows from representations (|4.13() , (|4.15|) 
and H4.17JI . We divide the left- and right-hand sides by t, and pass to the limits in 
probability as t increases to infinity. 
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Specifically, (|4.19jl and l|4.20jl follow from Ij4.13jl and from the fact that the 
martingales rriij(t) vanish with probability 1 as t — > oo. The last is the consequence 
of vanishing mn^i) and mu i+1 {t) with probability 1 as t — > oo. (The right-hand 
side of 14.14(1 is divided by t, and then we pass to the limit with probability 1 
as t — » oo.) Relations l(4.21[l and (|4.22(1 follow similarly from the corresponding 
relations l|4.15|l and (|4.17|l and vanishing the martingales m n fi(t) and rn n ^\{t) given 
by (|4.16l) and l|4.18|l . We also use the Lebesgue dominated convergence theorem to 
replace P- lim^^ by lim^oo E in the places where it is required. □ 

The standard particular cases from Theorem 14.11 are the cases where A(t) is 
Poisson, D(t) is Poisson, and both A(t) and D(t) are Poisson all together. We do 
not consider the first two cases, and the only last of these three cases is formulated 
and proved below. In this case the representation is deduced explicitly and has the 
form of limiting distributions rather than frequencies. 

Corollary 4.2. In the case where A(t) is Poisson with rate X, and D(t) is Poisson 
with rate A denote Py = lim^oo ¥{Qi(t) = i, Q2{t) = j}. We have the following. 
In the boundary cases i = for j = 0, 1 we have: 

(4.23) ///',., A/1, , • jAI'::,,. 

In the cases i = 1, 2, . . . , n — 1, j = 0, 1 we have: 

(4.24) ifiPij - (i + l)nPi+i,j = XPi-x,j - m, 3 + (1 - j)A/Vi,i - jPi,v 
In the case i — n and j — we have: 

(4.25) n/*P n , = AP„_i, - AP„, + AP n>1 . 
In the case i = n and j = 1 we have: 

(4.26) nf iP nA = AP n _i,i + AP„, . 

Proof. In the case where A(t) and D(t) are Poisson, the proof of this corollary is 
based on the following two relations: 

i r l 

lim -E / I{Qi(s-) = i, Q 2 (s-) = j}dA(s) = XP itj) 

(4.27) t_>0 ° ^ "'o 

i = 0,1, ...,n; j = 0,l, 
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and 

i r* 

lim -E / I{Qi(«-) = i,Q 2 (s-) = j}dD(s) = AP h3 , 
(4.28) 1 

i = 0, 1, ...,n; j = 0,1. 

Then inserting l|4.27j) and l|4.28j) into Ij4.19jl - I|4.22jl leads to the statement of the 
corollary. Thus there is only required to prove (|4.27|) and (|4.28|) . 

Prove (|4.27() . By semimartingale decomposition of Poisson process we have 
A(t) = Xt + Mji(t). Substituting this representation for (|4.27|l we obtain: 



(4.29) 



1 f* 

lim -E / I{Q 1 (s-)=i,Q 2 (s-)=j}dA(s) 

t-KX) t Jq 

a r* 

= lim - / F{Q 1 (s-)=i,Q 2 (s-)=j}ds 

t — oo t J Q 

1 /■* 

lim -E / I{Qi(s-) = i,Q 2 (s-)=j}dM A (s) 



t- 

i = 0, l,...,n; j = 0,1. 

The first term of the right-hand side is equal to XPi.j , and the second term is equal 
to 0. Therefore l|4.27|l is proved. 

The proof of (|4.28|) is analogous. The corollary is proved. □ 

5. Performance analysis and algorithms 

One of the most important performance characteristics is the loss proportion per 
arriving customer. Specifically, denoting the loss proportion by /, we write 

(5.1) P-hm^^^ l 

v ' 3 P-hmt^oo a{t) X w 

The parameters A, /i and A are assumed to be given. So, the problem is to find the 
appropriate number of servers n such that / < a. More accurately, the problem is 
to minimize n subject to / < a. 

In general case, we have no explicit representation for the processes, and then 
simulation techniques are used. There are three equivalent relations obtained for 
P-lim^oo^) in Section (E3J), (EH and (pHTIjl . 

The simplest of them looks (|3.1U|) . According to l|3.10fl the appropriate limit 
P- limt^+oo (73 (t) can be obtained straightforwardly. Immediately before each arrival 
(one occurring say at instant s) we only check the event {Qi(s— ) = n n Q2(s— ) = 
1}. Such type of simulation need not require any theory. However, in the case of 
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small a this simulation is not effective. As the large number of servers is chosen, 
and the loss proportion should be small, resulting in erroneous conclusion. 

The other relation is given by i|3.9|) . For a small a it is also based on the rare 
events, and its application is impossible. 

In this general case only (|3.8(1 is available. It is based on the average number of 
busy servers lim^oo t -1 E L Qi(s)ds which can be estimated by simulation. More 
specifically, the simulation is based on application of Theorem 14.11 helping us to 
obtain the frequencies 

lim - [ P{Qi(s) = i}ds = lim - / P{Qi(s) = i,Q 2 {s) = 0}ds 

i 

(5 - 2) + lim - / P{Qi(«) = i, Q 2 (s) = l}ds, 

t^OO t Jq 

i = 0, 1, . . . , n. 

The first and second terms of the right-hand size are obtained from Theorem 14. II 
They are expressed via the frequencies 

1 f* 

lim -E / I{Qi(«-) = i, Q 2 (s-) = j}dA(s) 

1 I* 

lim -E / I{Qi(s-) = i, Q 2 (s~) = j}dD(s) 

t^CO t Jq 

i = 0, 1, . . . ,n; j = 0,1, 

which in turn are easily evaluated by simulation. 

It is worth noting that direct simulating of lim^oo i _1 f Q P{Qi(s) = i, Q 2 {s) = 
j}ds (i = 0,1,..., n, j — 0,1) and consequently lim t ^ 00 t _1 E L Qi(s)ds is not 
available. (See Section 9 of for the detailed discussion of this question.) 

Thus, for the simulation purposes we use (|3.8() and the relations of Theorem 14. II 

Recall that the problem is to minimize m subject to / < a. The upper bound for 
n is unknown, and the special search procedure is necessary. The relevant search 
procedure is known due to Rubalskii |21j . who proposed the search algorithm for 
minimization of a unimodal function on an unbounded set. The optimal algorithm 
is an extension of the standard Fibonacci procedure. 

Thus, solution of the problem is based on two main steps: 

• Simulation 

• Search step 
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These procedures are repeated until the optimal solution is not found. 

In the particular case where both A(t) and D(t) are Poisson, for stationary 
probabilities Pij, i = 0, 1, . . . , n; j = 0, 1, we have the system of algebraic equations. 
Furthermore, the upper bound for the value n can be evaluated as well. In order 
to evaluate the upper bound for n, one can imagine that arrival rate is A + A, that 
is retrials occur all the time continuously and the retrial space is always occupied. 
In this case the loss probability is greater then that original, and consequently the 
corresponding value n is greater than that required. This upper bound n can be 
easily calculated by the Erlang loss formula. Then, having the upper and lower 
bounds we obtain the optimal number of servers easier. 

The similar approach can be used and for non-Markovian system. In the most 
practical cases where input stream is recurrent the upper bound for n can be easily 
evaluated by known formulae for the loss systems. 

However in general we have to use one or other special algorithm, say the method 
of |21j . in order to find optimal value of n. 

6. Numerical examples 

In this section we provide numerical solution for two systems. The first system 
is the Markovian system with A = 10, A = 2, /i = 1 and a = 0.0001 In the 
second example we assume that the increments of the processes A(t) and D(t) are 
deterministic with the same parameters. 

6.1. Example 1. Assuming that the lower bound for m is equal to 1, let us find 
the upper bound. For this purpose we consider the M/M/n/0 queueing system 
with input rate 12=A+A and /j, = 1. According to Erlang loss formula, for the loss 
probability we have 



Z^=o •, 

where p = A//i=12. Assuming that the loss probability is equal to a = 0.0001 let 
us find the value n. 

For n=26 the loss probability is 0.000174, and for n—27 the loss probability is 
0.000078 The first of these probabilities is greater than 0.0001, and the second one 
is smaller than 0.0001 Therefore, the upper bound for n is 27. 
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Let us now find n optimal. The solution of this problem can be achieved as 
follows. We find a new value of n=l+(27-l)/2=14. Now we solve the system of 
equations l|4.23|l - (|4.26|l given by Corollary 14.21 We also use (|5.1|l and l|3.13(l . For 
loss probability we obtain 0.064841 Therefore in the next step we find a new value 
n by 14+(27-14)/2=20.5 taking then the integer part. With new value n — 20 we 
solve the system again, and for the loss probability we obtain 0.001226 This value is 
greater than a, and therefore the new value of n is the integer part of 20+(27-20)/2 
equal to 23. With this value of n for the loss probability we have 0.000110 Thus 
the last step is anticipated with n=24. (Formally we must analyze first the value 
n=25 and only then accept n=24.) For this value n the loss probability is 0.000045 
Thus the problem is solved, and the value n=24 is the optimal number of servers 
decreasing the loss probability to the required value. 

6.2. Example 2. Assume that all parameters are the same as in Example 1, but 
increments of the processes A(t) and D(t) are deterministic. In this case the scheme 
of calculation includes simulation as it is explained in Sectional For this concrete 
system the algorithm can be simplified. Specifically, as in Example 1 we first 
evaluate the upper bound for n. For our approximation of the upper bound of 
n we consider the D/M/n/0 queueing system, the interarrival times of which all 
are equal to 1/(A + A), i.e. in our case 1/12. Notice, that under the assumption 
that customers constantly arrive to the main queue by A(t) and D(t) we do not 
have the model of D/M/n/0 queue exactly. Although input stream is based on 
two sources with deterministic interarrival times, the structure of the overall input 
stream is complicated. The arguments for such approximating by the D/M/n/0 
queueing system are heuristic. Nevertheless, such approximation is available. This 
is supported by computations as well. 

The loss probability for GI/M/n/0 loss systems is well-known (e.g. [B], JH], |20| . 



as well as 0): 




l 



(6.1) 



where 
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and R(x) is the probability distribution function of an interarrival time. In our case 
interarrival times are deterministic of length 1/12 and /j, = 1. Therefore, in (|6.1|1 

Tj =C-->/12 . 

According to our calculations the upper bound for n is 22. The value of the loss 
probability is 0.000036 Let us note that for n=21 the value of the loss probability 
is 0.000137, i.e. it is slightly greater than a=0.0001 

In the next step we choose then the integer part of l+(22-l)/2, and the value n is 
11. Simulating with this value of n yields the value of the loss probability 0.125444 
Therefore in the next step we choose the integer part of ll+(22-ll)/2, i.e. n=16. 
With this value of n by simulating we have the value of loss probability 0.002784 In 
the next step n=16+(22-16)/2=19. With this value n we correspondingly obtain 
the loss probability 0.000060 In the next step we choose the integer part of 19- 
(19-16)/2 for n, i.e. n=17. The loss probability in this case is 0.000786 There 
is one value n=18 which must be checked. For this value the loss probability is 
0.000214 Thus our conclusion is n=17. Recall that the loss probability in this case 
is 0.000060 

7. Concluding remarks 

In this paper we studied multiserver retrial queueing system having only one 
space in orbit. For this queueing system we derived analytic results permitting 
us to provide performance analysis. Specifically we solved the problem to find 
the possibly minimal number of servers decreasing the loss probability to given 
small value. In general case the algorithm of solution requires effective simulation, 
meeting the challenge of rare events problem. Numerical examples were provided for 
two cases. In the first case all processes were assumed to be Poisson, and the system 
Markovian. For Markovian system the calculations were based on analytic results. 
In the second case the processes A{t) and D(t) had deterministic increments. The 
analysis in this case was based on effective simulation. Comparison of these results 
showed that in the case of Markovian system the required number of servers, making 
the loss proportions smaller than the given value, is relatively greater than that in 
the system having the processes A(t) and D(t) with deterministic increments. 
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